<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta name="generator" content="MATLAB R2016a"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><title>Steady-state flow recirculation zone around a well doublet</title><style type="text/css">
* {margin: 0; padding: 0;}
body {text-align: start; line-height: 17.2339992523193px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Consolas, Inconsolata, Menlo, monospace; font-style: normal; font-size: 14px; font-weight: normal; text-decoration: none; white-space: normal; }
h1, h2 {font-weight: normal;}
.content { padding: 30px; }

.S0 { margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S1 { line-height: 26.3999996185303px; min-height: 24px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-size: 22px; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 15px; margin-right: 10px;  }
.S2 { min-height: 0px; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S3 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S4 { font-family: Helvetica, Arial, sans-serif; margin-left: 0px; margin-top: 10px; margin-bottom: 20px; margin-right: 0px;  }
.S5 { text-align: left; line-height: 21px; white-space: pre-wrap; white-space: pre-wrap; margin-left: 56px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S6 { line-height: 20.576000213623px; min-height: 20px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-size: 16px; font-weight: bold; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 9px; margin-right: 10px;  }
.S7 { margin-left: 3px; margin-top: 10px; margin-bottom: 10px; margin-right: 3px;  }
.S8 { line-height: 15.5926666259766px; min-height: 18px; white-space: nowrap; font-size: 12.6666669845581px; white-space: nowrap; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S9 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(34, 139, 34); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S10 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S11 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 10px; margin-bottom: 9px; margin-right: 10px;  }
.S12 { min-height: 0px; font-family: monospace; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S13 { min-height: 0px; color: rgb(0, 95, 206); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S14 { font-family: monospace; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S15 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S16 { margin-left: 3px; margin-top: 10px; margin-bottom: 4px; margin-right: 3px;  }
.S17 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(160, 32, 240); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S18 { color: rgb(64, 64, 64); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S19 { min-height: 0px; font-style: italic; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S20 { text-align: center; line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S21 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S22 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }

.LineNodeBlock {margin: 10px 0 10px 0;}
.LineNodeBlock+.paragraphNode {margin-top: 10px;}
.lineNode {padding-left: 10px; background-color: #F7F7F7; border-left: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.inlineWrapper:first-child .lineNode,.inlineWrapper.outputs+.inlineWrapper .lineNode {padding-top: 5px; border-top: 1px solid #E9E9E9;}
.inlineWrapper:last-child .lineNode,.inlineWrapper.outputs .lineNode {padding-bottom: 5px; border-bottom: 1px solid #E9E9E9;}
.lineNode .textBox {white-space: pre;}
.outputGroup {margin: 2px 0 2px 0; padding: 2px 2px 2px 4px;}
.outputRegion {}
.outputParagraph {color: rgba(64, 64, 64, 1); padding: 10px 0 6px 17px; background: white; overflow-x: hidden;}
.inlineWrapper:last-child .outputParagraph {border-bottom-left-radius: 4px; border-bottom-right-radius: 4px;}
.outputParagraph:empty {margin: 0;}
.inlineElement .symbolicElement {}
.embeddedOutputsSymbolicElement .MathEquation {margin-top: 5px; margin-bottom: 5px;}
.embeddedOutputsSymbolicElement .MathEquation.displaySymbolicElement {margin-left: 15px;}
.embeddedOutputsSymbolicElement .MathEquation.inlineSymbolicElement {margin-left: 5px;}
.embeddedOutputsSymbolicElement {overflow-x: auto; overflow-y: hidden;}
.embeddedOutputsSymbolicElement { overflow: initial !important;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre; word-wrap: initial; min-height: 18px; max-height: 250px; overflow: auto;}
.textElement {padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement { max-height: none !important; overflow: initial !important;}
.veSpecifier {}
.veContainer {}
.veSpecifierBox {height: 400px; width: 500px;}
.veSpecifier .veTable {padding-top: 3px; padding-bottom: 4px;}
.veSpecifierBox .veSpecifier .veContainer {position: relative; width: 100%; height: 370px;}
.veSpecifierBox .dijitDialogPaneContent {width: 97% !important; height: 88% !important;}
.veSpecifier .veTable .rowHeadersWrapper {padding-bottom: 0;}
.veSpecifier .veTable .scroller .variableEditorRenderers {padding-right: 3px; -webkit-user-select: none; -moz-user-select: none; -ms-user-select: none;}
.veVariableValueSummary {display: inline-block;}
.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {visibility: hidden; z-index: 0;}
.veMetaSummary {font-style: italic;}
.veSpecifier .veTable .scroller {overflow: hidden;}
.veSpecifier .veTable:hover .scroller {overflow: auto;}
.veSpecifier .veVariableName,.veSpecifier .veDimensionFont {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.veSpecifier .veVariableName {padding-top: 3px;}
.veSpecifier .veDimensionFont {font-style: italic; color: #9A9A9A;}
.veSpecifier .scroller::-webkit-scrollbar-track {background-color: white;}
.veSpecifier .scroller::-webkit-scrollbar-corner {background-color: white;}
.veSpecifier .veTable .topRowHeaderWrapper {border: none; background-color: #F8F9FA;}
.veSpecifier .mw_type_ListBase.showCellBorders,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper,.veSpecifier .veTable .verticalScrollSpacer,.veSpecifier .veTable .horizontalScrollSpacer {border: none;}
.veSpecifier .veTable .dataScrollerNode {border: 1px solid #BFBFBF;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode,.veSpecifier .veTable .dataBody {font-family: Arial; font-size: 13px;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode {color: #7F7F7F;}
.veSpecifier .veTable .dataBody {color: #000000;}
.veSpecifier .veTable .columnHeaderNode .cell .drag,.veSpecifier .veTable .columnHeaderNode .cell .header,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {background-color: #F8F9FA;}
.veSpecifier .veTable .columnHeaderNode .cell .dragBorder {border-right: 1px solid #F8F9FA;}
.veSpecifier .veTable .rowHeaderNode .cell {padding-top: 3px; text-align: center; border-bottom: 1px solid #F8F9FA;}
.veSpecifier .veTable .dataBody .cell .plainText {text-align: right;}
.veSpecifier .veTable .dataBody .row .cell {border-bottom: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.embeddedOutputsVariableElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre-wrap; word-wrap: break-word; min-height: 18px; max-height: 250px; overflow: auto;}
.variableElement {}
.embeddedOutputsVariableElement.inlineElement {}
.inlineElement .variableElement {}
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement {min-height: 18px; box-sizing: border-box; font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.matrixElement .variableValue {white-space: pre; display: inline-block; vertical-align: top; overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {padding-top: 0 !important; max-height: 100px;}
.embeddedOutputsMatrixElement.inlineElement .veVariableName {padding-top: 0;}
.inlineElement .matrixElement {max-height: 300px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer {white-space: nowrap; padding-top: 10px; margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide {display: none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis {margin-bottom: -3px;}
.matrixElement { max-height: none !important; }
.dijitTooltipDialog .dijitTooltipContainer .dijitTooltipContents .alertPlugin-alertMessage {min-width: 12px; max-width: 400px; max-height: 300px; overflow: auto;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar {width: 11px; height: 11px;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-track {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-corner {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {border-radius: 7px; background-color: rgba(0, 0, 0, 0.1); border: 2px solid rgba(0, 0, 0, 0); background-clip: padding-box;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0.1);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.alertPlugin-alertLine {position: absolute; display: initial; width: 40px; text-align: right; cursor: text;}
.alertPlugin-onTextLine {visibility: hidden;}
.alertPlugin-hasTooltip .alertPlugin-warningImg,.alertPlugin-hasTooltip .alertPlugin-errorImg {cursor: pointer;}
.alertPlugin-isStale {-webkit-filter: opacity(0.4) grayscale(80%); filter: opacity(0.4) grayscale(80%);}
.alertPlugin-alertLine .alertPlugin-errorElement {display: inline-block; margin-right: 4px;}
.alertPlugin-errorImg {position: relative; display: inline-block; background-color: rgb(204, 55, 41); width: 14px; height: 14px; margin-top: 1px; -webkit-border-radius: 7px;}
.alertPlugin-alertLine .alertPlugin-warningElement {display: inline-block; margin-right: 3px;}
.alertPlugin-warningImg {position: relative; display: inline-block; width: 16px; height: 15px; overflow: hidden;}
.diagnosticMessage-wrapper {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {color: rgb(255,100,0); text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {color: rgb(230,0,0); text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart {white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {white-space: pre;}
.embeddedOutputsWarningElement{min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsErrorElement {min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.alertPlugin-warningImg {background-image: url()}
.alertPlugin-errorImg {background-image: url()}
</style></head><body><div class = "content"><div class = 'SectionBlock containment'><h1 class = "S1"><span class = "S2">Steady-state flow recirculation zone around a well doublet</span></h1><p class = "S3"><span class = "S2">Flow recirculation zones around one or many pairs of injection-extraction wells occur in many practical groundwater flow, reservoir engineering and geothermal engineering problems. Typical application examples are capture zone delineation, forced-gradient tracer tests for aquifer characterization, in-situ remediation of contaminated groundwater, thermogeological assessment of open-loop well-doublet schemes, seasonal heat storage and recovery, and cold fluid reinjection into produced geothermal reservoirs. </span></p><p class = "S3"><span class = "S2">The purpose of this modelling exercise is to simulate steady-state flow patterns around a pair of injection and production wells which are supposed to inject back pumped water with the same flow rate. </span></p><p class = "S3"><span class = "S2">Problem dimensions and aquifer characteristics used in this tutorial are patterned following the mean characteristics of the carbonates Dogger geothermal reservoir layer which is a mildly 1500 meters deep hot aquifer in the Paris basin. It represents the only operational low-enthalpy geothermal system in France since the 1980's. </span></p><p class = "S3"><span class = "S2">We will go through the following steps to build a simulation script modelling the steady-state flow patterns around one doublet in this geothermal reservoir:</span></p><ol class = "S4"><li class = "S5"><span class = "S0">Setup the computational grid, wells positions, flow rates, and fluid properties.</span></li><li class = "S5"><span class = "S0">Perform a steady-state flow simulation around the well pair. </span></li><li class = "S5"><span class = "S0">Visualize the overpressure distribution in the recirculating zone. </span></li><li class = "S5"><span class = "S0">Compute and visualize flow streamlines in this zone. </span></li><li class = "S5"><span class = "S0">Compare numerically computed streamlines with analytical solutions derived from the complex potential theory for two-dimensional homogeneous and isotropic Dupuit-Forcheimer aquifers. </span></li></ol></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Computational grid, wells, and fluid properties setup </span></h2><p class = "S3"><span class = "S2">Let's first specify the simple aquifer geometry which is a one layer model of 6000m x 6000m x 25m along each space dimension.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Domain size along X, Y &amp; Z directions</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Dx = 6000; Dy = 6000; Dz = 25;</span></p></div></div><p class = "S11"><span class = "S2">Next, we set the number of grid cells along each direction. Herein we choose a uniform grid spacing equal to 25m along X and Y directions in the horizontal plane. Next, let </span><span class = "S12">N</span><span class = "S2"> to equal the number of grid cells (or elements) in the model.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Nx = 240;  Ny = 240;  Nz = 1;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">N  = Nx*Ny*Nz;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.Nx = Nx; Grid.Ny = Ny; Grid.Nz = Nz; Grid.N = N;</span></p></div></div><p class = "S11"><span class = "S2">To fullfill the need of internal compuattional tasks we add the uniform grid spacings, </span><span class = "S12">hx</span><span class = "S2">, </span><span class = "S12">hy</span><span class = "S2">, and </span><span class = "S12">hz</span><span class = "S2">, along the respective directions to the </span><span class = "S12">Grid</span><span class = "S2"> structure</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hx = (Dx/Nx); Grid.hy = (Dy/Ny); Grid.hz = (Dz/Nz);</span></p></div></div><p class = "S11"><span class = "S2">Additionally, we need to setup other few variables for the </span><span class = "S12">Grid</span><span class = "S2"> structure such as, </span><span class = "S12">V</span><span class = "S2"> (the uniform cells volume), and </span><span class = "S12">K</span><span class = "S2"> (the tintrinsic permeability tensor). Notice that the permeability data is a 3D shape array whose size is 3*Nx*Ny*Nz. The MATLAB command </span><a href = 'https://fr.mathworks.com/help/matlab/ref/ones.html'><span class = "S14">ones</span></a><span class = "S12">(3,Nx,Ny,Nz)</span><span class = "S2"> constructs such an array with a unit values for all cells, and thus by multiplying it with a constant unit Darcy we obtain an isotropic and homogeneous distribution of the permeability data over the computational domain. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">darcy = 9.869233e-13;         </span><span class = "S9">% factor for conversion from Darcy to SI units </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.K = darcy*ones(3,Nx,Ny,Nz);</span></p></div></div><p class = "S11"><span class = "S2">Next, we construct a cartesian grid </span><span class = "S12">G</span><span class = "S2"> data structure from the domain extents and spacings. This will help to do the plotting commands more efficiently.   </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);</span></p></div></div><p class = "S11"><span class = "S2">Before setting the flow rates we define some useful constants to convert between time units</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">second = 1; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">minute = 60*second; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">hour   = 60*minute; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">day    = 24*hour; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">month  = 30*day;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">year   = 365*day;</span></p></div></div><p class = "S11"><span class = "S2">We set the injection and production flow rates equal to +150</span><span style="vertical-align:-3"><img src="" width="36" height="24" /></span><span class = "S2"> and -150</span><span style="vertical-align:-3"><img src="" width="36" height="24" /></span><span class = "S2">, respectively. Notice that (by convention) injection flow rates are positive and production ones are negrative. </span><span class = "S12">Qw</span><span class = "S2"> is the well rates array. It is null everywhere except at cell numbers where the wells are located. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">wells_pos = [91,150];              </span><span class = "S9">% cell no's of the injection and production wells respectively</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">Q  = 150/hour;                     </span><span class = "S9">% conversion from m^3/h to m^3/s</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">Qw = zeros(N,1);                   </span><span class = "S9">% initialization</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">Qw(wells_pos) = Q*[+1,-1];         </span><span class = "S9">% injection/production wells flow rates</span></p></div></div><p class = "S11"><span class = "S2">Finally, we setup the fluid properties by initializing a </span><span class = "S12">Fluid</span><span class = "S2"> object named 'water' with user-defined constant density and viscosity properties equal to </span><span style="vertical-align:-3"><img src="" width="67" height="24" /></span><span class = "S2"> and </span><span style="vertical-align:-3"><img src="" width="56.5" height="21" /></span><span class = "S2"> respectively. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% fluid properties units: density [Kg/m^3], viscosity [Pa.s]</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">water = Fluid(</span><span class = "S17">'water'</span><span class = "S10">,[1000, 1e-3]);</span></p></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Steady-state flow simulation</span></h2><p class = "S3"><span class = "S2">To simulate steady-state flow we simply call the </span><span class = "S12">Pressure</span><span class = "S2"> solver. In this flow regime this function needs only three input arguments: the computational grid, </span><span class = "S12">G3</span><span class = "S2">, the fluid, </span><span class = "S12">water</span><span class = "S2">, and the array of flow rates, </span><span class = "S12">Qw</span><span class = "S2">. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Pascal = 1;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">bar    = 1e5*Pascal;               </span><span class = "S9">% conversion from bar to Pascal units</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% call flow solver for single-phase flow</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S10">tic; [P,V] = Pressure(Grid,water,Qw); toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="936" style="width: 936px;"><div class="textElement">Elapsed time is 0.747978 seconds.</div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Overpressure distribution in the recirculation zone</span></h2><p class = "S3"><span class = "S2">The overpressure distribution arround the wells doublet is plotted using the </span><span class = "S12">PlotContourData</span><span class = "S2"> function. Here, we call </span><a href = 'https://fr.mathworks.com/help/matlab/ref/figure.html'><span class = "S14">figure</span></a><span class = "S2"> command to open a new figure window so that all next plotting commands will be directed into it. The plot is presented with contour levels between the min and max levels of overpressure around </span><span class = "S19">-25 bar</span><span class = "S2"> and 25</span><span class = "S19"> bar</span><span class = "S2"> respectively with 1.25 bar interval. We use the </span><a href = 'https://fr.mathworks.com/help/matlab/ref/colormap.html'><span class = "S14">jet</span></a><span class = "S2"> </span><a href = 'https://fr.mathworks.com/help/matlab/ref/colormap.html'><span class = "S14">colormap</span></a><span class = "S2"> with 41 levels of colors equating thus the number of intervals between contour lines, thus producing a plot in which each region between two succesive contour lines is filled with a unique color. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot overpressure in the domain</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">opts.fill = </span><span class = "S17">'off'</span><span class = "S10">;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">PlotContourData(G, P(:)/bar, opts,</span><span class = "S17">'LevelList'</span><span class = "S10">,-25:1.25:25);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(41));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">axis </span><span class = "S17">tight equal</span><span class = "S10">; axis([0,6000,0,3000]); colorbar; </span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S15">title(</span><span class = "S17">'Overpressure contour lines [bar]'</span><span class = "S15">); xlabel(</span><span class = "S17">'X[m]'</span><span class = "S15">); ylabel(</span><span class = "S17">'Y[m]'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 936px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S11"><span class = "S2">As expected the flow pattern is perfectly symmetrical around the median y-axis of the rectangular domain and the contour lines are orthogonal to no-flow boundaries according to the classical groundwater flownet theory. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Streamlines in the recirculation zone</span></h2><p class = "S3"><span class = "S2">Now, we will calculate and plot the streamlines (lines tangent to velocities) starting from an equal distribution of points along the median y axis. We set massless particle tracking directions forwards towards the producer and backwards towards the injector. Next, we plot the two different sets of forward and backward streamlines with different colors. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">h = [Dx/Nx, Dy/Ny];                        </span><span class = "S9">% uniform grid sizes along X &amp; Y</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">sx = 3000*ones(1,7);                       </span><span class = "S9">% streamlines starting x positions</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">sy = linspace(0.5*h(2),(Ny-0.5)*h(2)/2,7); </span><span class = "S9">% streamlines starting y positions</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">hold </span><span class = "S17">on</span><span class = "S15">;                                   </span><span class = "S9">% keep drawing to last plot</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">hp = PlotStreamline(V,h,sx,sy,true);       </span><span class = "S9">% plot streamlines in forward direction ...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">set(hp,</span><span class = "S17">'Color'</span><span class = "S15">,</span><span class = "S17">'b'</span><span class = "S15">,</span><span class = "S17">'LineWidth'</span><span class = "S15">,0.8);       </span><span class = "S9">%    with blue color </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">hi = PlotStreamline(V,h,sx,sy,false);      </span><span class = "S9">% plot streamlines in backward direction ... </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">set(hi,</span><span class = "S17">'Color'</span><span class = "S15">,</span><span class = "S17">'r'</span><span class = "S15">,</span><span class = "S17">'LineWidth'</span><span class = "S15">,0.8);       </span><span class = "S9">%    with red color</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S15">title(</span><span class = "S17">'Overpressure contour lines and streamlines'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 936px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S11"><span class = "S2">Here again, the obtained streamlines are orthogonal to pressure (and therefore discharge potential) contour lines as expected from the flownet theory for stedy-state two-dimensional groundwater flow in homogeneous and isotropic aquifers. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment active'><h2 class = "S6"><span class = "S2">Comparison to analytical solution</span></h2><p class = "S3"><span class = "S2">When considering groundwater flow in homogeneous and isotropic confined aquifers, Bear (1979) and Strack (1989) have used the complex potential theory to analytically derive the expressions for its discharge potential, </span><span style="vertical-align:-3"><img src="" width="14.5" height="18" /></span><span class = "S2">, and the streamfunction, </span><span style="vertical-align:-3"><img src="" width="14.5" height="18" /></span><span class = "S2">, components assuming the Dupuit-Forcheimer conditions. The complex potential is given as</span></p><p class = "S20"><span style="vertical-align:-3"><img src="" width="132.5" height="20" /></span></p><p class = "S3"><span class = "S2">where </span><span style="vertical-align:-3"><img src="" width="62" height="18" /></span><span class = "S2"> is the coordinate in the complex plane. </span></p><p class = "S3"><span class = "S2">In the recirculation zone created by a pair of of injection and extraction wells in absence of regional flow these functions are given as </span></p><p class = "S20"><span style="vertical-align:-13"><img src="" width="156" height="40" /></span></p><p class = "S20"><span style="vertical-align:-9"><img src="" width="253.5" height="26" /></span></p><p class = "S3"><span class = "S2">where </span><span style="vertical-align:-3"><img src="" width="12.5" height="18" /></span><span class = "S2">is the half-distance between the wells.</span></p><p class = "S3"><span class = "S2">The equipotential lines are concentric circles centered at the injection and production wells. Their analytic expressions are given by Strack (1989) as</span></p><p class = "S3"><span style="vertical-align:-22"><img src="" width="215.5" height="45" /></span><span class = "S2">   if  </span><span style="vertical-align:-3"><img src="" width="41" height="18" /></span></p><p class = "S3"><span style="vertical-align:-3"><img src="" width="36" height="18" /></span><span class = "S2">                                                  if </span><span style="vertical-align:-3"><img src="" width="41" height="18" /></span></p><p class = "S3"><span class = "S2">The streamlines have also circular shapes and their center are located in the median vertical axis between the well pair. They're analytically expressed as </span></p><p class = "S3"><span style="vertical-align:-22"><img src="" width="201" height="45" /></span><span class = "S2">   if </span><span style="vertical-align:-9"><img src="" width="64.5" height="26" /></span><span class = "S2"> </span></p><p class = "S3"><span style="vertical-align:-3"><img src="" width="38.5" height="18" /></span><span class = "S2">                                             if </span><span style="vertical-align:-9"><img src="" width="64.5" height="26" /></span></p><p class = "S3"><span class = "S2">It is easy to show that coordinates of the streamline, </span><span style="vertical-align:-3"><img src="" width="14.5" height="18" /></span><span class = "S2">, center are </span><span style="vertical-align:-3"><img src="" width="118.5" height="20" /></span><span class = "S2">and that its radius is </span><span style="vertical-align:-3"><img src="" width="110.5" height="20" /></span><span class = "S2">. The height of streamline arc is </span><span style="vertical-align:-9"><img src="" width="120.5" height="28" /></span><span class = "S2">(Luo and Kitanidis, 2004).</span></p><p class = "S3"><span class = "S2">Now, let's calculate the </span><span style="vertical-align:-3"><img src="" width="14.5" height="18" /></span><span class = "S2"> levels of numerically calculated streamlines from their prescribed heights. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">d   = 1500/2;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">psi = Q*((pi/2)-atan(sy(2:end)/d))/pi;</span></p></div></div><p class = "S11"><span class = "S2">Now, we plot few points of the corresponding circles for each streamline except the first one corresponding to the straight line connecting the two wells. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">rc = d./abs(sin(2*pi*psi/Q));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">yc = d*cot(2*pi*psi/Q);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S21">for </span><span class = "S10">i=1:numel(rc)</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    xp = rc(i)*cos(-pi:pi/21:pi)+Dx/2; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    yp = rc(i)*sin(-pi:pi/21:pi)+yc(i); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S15">    plot(xp,yp,</span><span class = "S17">'bo'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S22">end</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S15">title(</span><span class = "S17">'Overpressure, analytical and numerical streamlines'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 936px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S11"><span class = "S2">As shown in the last figure the agreement between the numerically computed streamlines (plotted in solid lines) and the analytically calculated ones (represented by circles) is quite good. The quality of this agreement deteriorates for the longest streamlines connecting the wells. This is a boundary condition artifact since those streamlines become closer to the closed domain boundaries, unlike in the analytical derivations assuming an infinite domain extent. This error should not be attributed to numerical issues associated with the pressure solver or the streamlines integrator. This could be easily checked by decreasing the grid size in both directions showing to have negligeable impact on the numerical solution accuracy. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Conclusions</span></h2><p class = "S3"><span class = "S2">In this tutorial we learned how to make a simple script to design the recirculation zone around a well doublet in a geothermal reservoir to establish the initial steady-state flow patterns prior to a heat transport simulation. </span></p><p class = "S3"><span class = "S2">Numerically calculated and visualized flow streamlines delineating the recirculation zone are compared to analytic expressions derived from the complex potential theory. A good agreement is observed between these two solutions demonstrating the accuracy of the flow solver and the streamlines tracer. </span></p><p class = "S3"><span class = "S2">This tutorial could be easily adapted to jump-start other applications. For instance, in capture zone delineation problems the variable of interest is the volume of each recirculation zone around one doublet or a complex system of hydraulically connected wells. Otherwise, for in-situ remediation or tracer test problems we need to consider the residence or arrival times. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">References</span></h2><p class = "S3"><span class = "S2">Bear, J. (1979). </span><span class = "S19">Hydraulics of Groundwater</span><span class = "S2">. McGraw-Hill, New York. </span></p><p class = "S3"><span class = "S2">Luo, J., and P.K. Kitanidis (2004). Fluid residence times within a recirculation zone created by an extraction-injection well pair. </span><span class = "S19">J. Hydrology</span><span class = "S2">, 295, 149-162. </span></p><p class = "S3"><span class = "S2">Strack, O.D.L. (1989). </span><span class = "S19">Groundwater Mechanics</span><span class = "S2">, Prentice-Hall, Englewood Cliffs, NJ. </span></p></div></div>
<!-- 
##### SOURCE BEGIN #####
%% Steady-state flow recirculation zone around a well doublet
% Flow recirculation zones around one or many pairs of injection-extraction 
% wells occur in many practical groundwater flow, reservoir engineering and geothermal 
% engineering problems. Typical application examples are capture zone delineation, 
% forced-gradient tracer tests for aquifer characterization, in-situ remediation 
% of contaminated groundwater, thermogeological assessment of open-loop well-doublet 
% schemes, seasonal heat storage and recovery, and cold fluid reinjection into 
% produced geothermal reservoirs. 
% 
% The purpose of this modelling exercise is to simulate steady-state flow 
% patterns around a pair of injection and production wells which are supposed 
% to inject back pumped water with the same flow rate. 
% 
% Problem dimensions and aquifer characteristics used in this tutorial are 
% patterned following the mean characteristics of the carbonates Dogger geothermal 
% reservoir layer which is a mildly 1500 meters deep hot aquifer in the Paris 
% basin. It represents the only operational low-enthalpy geothermal system in 
% France since the 1980's. 
% 
% We will go through the following steps to build a simulation script modelling 
% the steady-state flow patterns around one doublet in this geothermal reservoir:
% 
% # Setup the computational grid, wells positions, flow rates, and fluid properties.
% # Perform a steady-state flow simulation around the well pair. 
% # Visualize the overpressure distribution in the recirculating zone. 
% # Compute and visualize flow streamlines in this zone. 
% # Compare numerically computed streamlines with analytical solutions derived 
% from the complex potential theory for two-dimensional homogeneous and isotropic 
% Dupuit-Forcheimer aquifers. 
%% Computational grid, wells, and fluid properties setup 
% Let's first specify the simple aquifer geometry which is a one layer model 
% of 6000m x 6000m x 25m along each space dimension.

% Domain size along X, Y & Z directions
Dx = 6000; Dy = 6000; Dz = 25;
%% 
% Next, we set the number of grid cells along each direction. Herein we 
% choose a uniform grid spacing equal to 25m along X and Y directions in the horizontal 
% plane. Next, let |N| to equal the number of grid cells (or elements) in the 
% model.

Nx = 240;  Ny = 240;  Nz = 1;
N  = Nx*Ny*Nz;
Grid.Nx = Nx; Grid.Ny = Ny; Grid.Nz = Nz; Grid.N = N;
%% 
% To fullfill the need of internal compuattional tasks we add the uniform 
% grid spacings, |hx|, |hy|, and |hz|, along the respective directions to the 
% |Grid| structure

Grid.hx = (Dx/Nx); Grid.hy = (Dy/Ny); Grid.hz = (Dz/Nz);
%% 
% Additionally, we need to setup other few variables for the |Grid| structure 
% such as, |V| (the uniform cells volume), and |K| (the tintrinsic permeability 
% tensor). Notice that the permeability data is a 3D shape array whose size is 
% 3*Nx*Ny*Nz. The MATLAB command |<https://fr.mathworks.com/help/matlab/ref/ones.html 
% ones>(3,Nx,Ny,Nz)| constructs such an array with a unit values for all cells, 
% and thus by multiplying it with a constant unit Darcy we obtain an isotropic 
% and homogeneous distribution of the permeability data over the computational 
% domain. 

darcy = 9.869233e-13;         % factor for conversion from Darcy to SI units 
Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);
Grid.K = darcy*ones(3,Nx,Ny,Nz);
%% 
% Next, we construct a cartesian grid |G| data structure from the domain 
% extents and spacings. This will help to do the plotting commands more efficiently.   

G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);
%% 
% Before setting the flow rates we define some useful constants to convert 
% between time units

second = 1; 
minute = 60*second; 
hour   = 60*minute; 
day    = 24*hour; 
month  = 30*day;
year   = 365*day;
%% 
% We set the injection and production flow rates equal to +150$m^3/h$ and 
% -150$m^3/h$, respectively. Notice that (by convention) injection flow rates 
% are positive and production ones are negrative. |Qw| is the well rates array. 
% It is null everywhere except at cell numbers where the wells are located. 

wells_pos = [91,150];              % cell no's of the injection and production wells respectively
Q  = 150/hour;                     % conversion from m^3/h to m^3/s
Qw = zeros(N,1);                   % initialization
Qw(wells_pos) = Q*[+1,-1];         % injection/production wells flow rates
%% 
% Finally, we setup the fluid properties by initializing a |Fluid| object 
% named 'water' with user-defined constant density and viscosity properties equal 
% to $10^3Kg/m^3$ and $10^{-3}Pa.s$ respectively. 

% fluid properties units: density [Kg/m^3], viscosity [Pa.s]
water = Fluid('water',[1000, 1e-3]);
%% Steady-state flow simulation
% To simulate steady-state flow we simply call the |Pressure| solver. In this 
% flow regime this function needs only three input arguments: the computational 
% grid, |G3|, the fluid, |water|, and the array of flow rates, |Qw|. 

Pascal = 1;
bar    = 1e5*Pascal;               % conversion from bar to Pascal units

% call flow solver for single-phase flow
tic; [P,V] = Pressure(Grid,water,Qw); toc;
%% Overpressure distribution in the recirculation zone
% The overpressure distribution arround the wells doublet is plotted using the 
% |PlotContourData| function. Here, we call |<https://fr.mathworks.com/help/matlab/ref/figure.html 
% figure>| command to open a new figure window so that all next plotting commands 
% will be directed into it. The plot is presented with contour levels between 
% the min and max levels of overpressure around _-25 bar_ and 25_ bar_ respectively 
% with 1.25 bar interval. We use the |<https://fr.mathworks.com/help/matlab/ref/colormap.html 
% jet>| |<https://fr.mathworks.com/help/matlab/ref/colormap.html colormap>| with 
% 41 levels of colors equating thus the number of intervals between contour lines, 
% thus producing a plot in which each region between two succesive contour lines 
% is filled with a unique color. 

figure; 
% plot overpressure in the domain
opts.fill = 'off';
PlotContourData(G, P(:)/bar, opts,'LevelList',-25:1.25:25);
colormap(jet(41));
axis tight equal; axis([0,6000,0,3000]); colorbar; 
title('Overpressure contour lines [bar]'); xlabel('X[m]'); ylabel('Y[m]');
%% 
% As expected the flow pattern is perfectly symmetrical around the median 
% y-axis of the rectangular domain and the contour lines are orthogonal to no-flow 
% boundaries according to the classical groundwater flownet theory. 
%% Streamlines in the recirculation zone
% Now, we will calculate and plot the streamlines (lines tangent to velocities) 
% starting from an equal distribution of points along the median y axis. We set 
% massless particle tracking directions forwards towards the producer and backwards 
% towards the injector. Next, we plot the two different sets of forward and backward 
% streamlines with different colors. 

h = [Dx/Nx, Dy/Ny];                        % uniform grid sizes along X & Y
sx = 3000*ones(1,7);                       % streamlines starting x positions
sy = linspace(0.5*h(2),(Ny-0.5)*h(2)/2,7); % streamlines starting y positions
hold on;                                   % keep drawing to last plot
hp = PlotStreamline(V,h,sx,sy,true);       % plot streamlines in forward direction ...
set(hp,'Color','b','LineWidth',0.8);       %    with blue color 
hi = PlotStreamline(V,h,sx,sy,false);      % plot streamlines in backward direction ... 
set(hi,'Color','r','LineWidth',0.8);       %    with red color
title('Overpressure contour lines and streamlines');
%% 
% Here again, the obtained streamlines are orthogonal to pressure (and therefore 
% discharge potential) contour lines as expected from the flownet theory for stedy-state 
% two-dimensional groundwater flow in homogeneous and isotropic aquifers. 
%% Comparison to analytical solution
% When considering groundwater flow in homogeneous and isotropic confined aquifers, 
% Bear (1979) and Strack (1989) have used the complex potential theory to analytically 
% derive the expressions for its discharge potential, $\Phi$, and the streamfunction, 
% $\Psi$, components assuming the Dupuit-Forcheimer conditions. The complex potential 
% is given as
% 
% $$\Omega(z) = \Phi(z) + i \Psi(z)$$
% 
% where $z=x+iy$ is the coordinate in the complex plane. 
% 
% In the recirculation zone created by a pair of of injection and extraction 
% wells in absence of regional flow these functions are given as 
% 
% $$\Phi(x,y) = \frac{Q}{4\pi} \ln [\frac{(x-d)^2+y^2}{(x+d)^2+y^2}]$$
% 
% $$\Psi(x,y) = \frac{Q}{2\pi}[\arctan(\frac{y}{x-d})-\arctan (\frac{y}{x+d})]$$
% 
% where $d$is the half-distance between the wells.
% 
% The equipotential lines are concentric circles centered at the injection 
% and production wells. Their analytic expressions are given by Strack (1989) 
% as
% 
% $[x+d\coth(\frac{2\pi\Phi}{Q})]^2 + y^2 = \frac{d^2}{\sinh^2(\frac{2\pi\Phi}{Q})}$   
% if  $\Phi \ne 0$
% 
% $x=0$                                                  if $\Phi = 0$
% 
% The streamlines have also circular shapes and their center are located 
% in the median vertical axis between the well pair. They're analytically expressed 
% as 
% 
% $x^2 + [y-d\cot(\frac{2\pi\Psi}{Q})]^2 = \frac{d^2}{sin^2(\frac{2\pi\Psi}{Q})}$   
% if $\Psi \ne 0, \pm\frac{Q}{2}$ 
% 
% $y=0$                                             if $\Psi = 0, \pm\frac{Q}{2}$
% 
% It is easy to show that coordinates of the streamline, $\Psi$, center are 
% $[0,d\cot(2\pi\Psi/Q)]$and that its radius is $d/|\sin(2\pi\Psi/Q)|$. The height 
% of streamline arc is $h_s = d\tan(\frac{\pi}{2}-\frac{\pi|\psi|}{Q})$(Luo and 
% Kitanidis, 2004).
% 
% Now, let's calculate the $\Psi$ levels of numerically calculated streamlines 
% from their prescribed heights. 

d   = 1500/2;
psi = Q*((pi/2)-atan(sy(2:end)/d))/pi;
%% 
% Now, we plot few points of the corresponding circles for each streamline 
% except the first one corresponding to the straight line connecting the two wells. 

rc = d./abs(sin(2*pi*psi/Q));
yc = d*cot(2*pi*psi/Q);
for i=1:numel(rc)
    xp = rc(i)*cos(-pi:pi/21:pi)+Dx/2; 
    yp = rc(i)*sin(-pi:pi/21:pi)+yc(i); 
    plot(xp,yp,'bo');
end
title('Overpressure, analytical and numerical streamlines');
%% 
% As shown in the last figure the agreement between the numerically computed 
% streamlines (plotted in solid lines) and the analytically calculated ones (represented 
% by circles) is quite good. The quality of this agreement deteriorates for the 
% longest streamlines connecting the wells. This is a boundary condition artifact 
% since those streamlines become closer to the closed domain boundaries, unlike 
% in the analytical derivations assuming an infinite domain extent. This error 
% should not be attributed to numerical issues associated with the pressure solver 
% or the streamlines integrator. This could be easily checked by decreasing the 
% grid size in both directions showing to have negligeable impact on the numerical 
% solution accuracy. 
%% Conclusions
% In this tutorial we learned how to make a simple script to design the recirculation 
% zone around a well doublet in a geothermal reservoir to establish the initial 
% steady-state flow patterns prior to a heat transport simulation. 
% 
% Numerically calculated and visualized flow streamlines delineating the 
% recirculation zone are compared to analytic expressions derived from the complex 
% potential theory. A good agreement is observed between these two solutions demonstrating 
% the accuracy of the flow solver and the streamlines tracer. 
% 
% This tutorial could be easily adapted to jump-start other applications. 
% For instance, in capture zone delineation problems the variable of interest 
% is the volume of each recirculation zone around one doublet or a complex system 
% of hydraulically connected wells. Otherwise, for in-situ remediation or tracer 
% test problems we need to consider the residence or arrival times. 
%% References
% Bear, J. (1979). _Hydraulics of Groundwater_. McGraw-Hill, New York. 
% 
% Luo, J., and P.K. Kitanidis (2004). Fluid residence times within a recirculation 
% zone created by an extraction-injection well pair. _J. Hydrology_, 295, 149-162. 
% 
% Strack, O.D.L. (1989). _Groundwater Mechanics_, Prentice-Hall, Englewood 
% Cliffs, NJ.
##### SOURCE END #####
--></body></html>